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Failure  Prediction  of  Underwater  Structures- 
Subdomain  Decomposition  and  Meshfree  Methods 


By  Ted  Belytschko 


Introduction: 

Meshfree  methods  for  fracture  have  been  extended  by  considering  more  general  classes  of  basis 
functions.  These  methods  do  not  require  any  elements,  and  boundaries  and  interfaces,  such  as  cracks, 
are  easy  to  propagate  in  the  model  since  it  does  not  involve  remeshing.  The  versatility  of  these  methods 
has  been  enhanced  by  developing  new  vector  level  set  methods.  In  these  methods,  the  crack  is 
represented  entirely  by  vector  functions  at  the  nodes  of  the  meshless  model.  Therefore,  there  is  no  need 
for  alternative  representations  of  the  crack  as  it  evolves.  This  extra  data  is  only  needed  at  nodes  in  a 
subdomain  around  the  crack.  The  method  has  been  applied  to  a  variety  of  crack  growth  problems  in  two 
dimensions  including  problems  involving  welds  and  fillets.  Comparisons  with  experiments  show  excellent 
agreement.  Extensions  of  the  method  to  three  dimensions  are  now  under  development.  Further  studies 
have  been  made  of  the  stability  of  particle  forms  of  meshless  methods.  It  has  been  shown  that  in  two 
dimensions  the  widely  used  Eulerian  kernels  distort  the  onset  of  material  instability  and  often  manifest 
instabilities  that  look  like  material  instabilities.  It  has  also  been  shown  that  Lagrangian  kernels  that  were 
proposed  earlier  in  this  grant  do  not  distort  the  onset  of  material  instabilities. 

Domain  decomposition  methods,  whereby  a  fine  scale  model  can  be  linked  with  a  coarse  scale  model,  are 
being  developed.  These  methods  are  useful  when  a  part  of  a  ship  where  severe  damage  is  expected  must 
be  modeled  by  a  refined  mesh  for  accuracy  and  linked  to  a  coarser  mesh  of  the  remainder  of  the  ship. 

The  method  that  has  been  developed  is  based  on  an  overlapping  domain  decomposition  method.  Part  of 
the  fine  scale  subdomain  overlaps  the  coarse  scale  subdomain.  In  each  time  step,  the  fine  scale  solution 
is  projected  on  the  coarse  scale  solution  so  that  spurious  reflections  on  the  interface  are  prevented. 
Methods  have  also  been  developed  in  which  a  smaller  time  step  can  be  used  in  the  fine  scale  subdomain. 
The  method  has  been  studied  in  model  one  and  two-dimensional  problems.  Results  show  that  spurious 
reflections  are  almost  totally  eliminated. 


1.  Stability  of  Meshless  Methods  based  on  Lagrangian  and  Eulerian  kernels 


A  major  difficult  of  meshless  methods  such  as  SPH  and  EFG  has  been  the  appearance  of  a  so-called 
tensile  instability.  We  have  shown  that  the  tensile  instability  is  even  more  pervasive  than  it  appears  at  first, 
for  it  drastically  changes  the  stable  material  domain.  Therefore,  the  prediction  of  cracks  or  shear  bands 
by  a  meshless  method  with  a  distorted  material  stability  domain  can  lead  to  highly  erroneous  results.  In 
the  following  we  show  the  analysis  which  has  led  us  to  that  conclusion.  We  also  propose  a  remedy,  the 
use  of  Lagrangian  kernels. 


An  isotropic  hyperelastic  material  model  is  used  the  study  the  domain  of  material  stability.  Its  potential  is 
given  by 


1  ,  i  r 

¥  =  ~^C\I\  +'^C2I2  ~ 


(1.1) 


where  c1  and  c2  are  constants,  /j ,  I2  and  /3  are  the  principal  invariants  of  the  right  Cauchy-Green 
deformation  tensor  C . 
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Based  on  a  linearized  stability  analysis  the  material  instability  in  the  continuum  is  given.  Fig.  1.1  shows  the 
stable  domain  of  continuum.  We  call  the  stable  domain  as  calculated  for  the  partial  differential  equation 
the  “exacf  stable  domain,  even  though  it  is  obtained  numerically. 


Fig.  1.1  :  The  stable  domain  of  continuum 


(a)  Lagrangian  kernel  (b)  Eulerian  kernel 

Fig.  1.2  :  Stable  domain  of  particle  methods  with  stress  point  integration 

We  studied  the  material  instability  for  particle  methods  with  either  Lagrangian  kernel  or  Eulerian  kernel. 
Fig.  1.2(a)  shows  that  Lagrangian  kernel  can  exactly  reproduce  the  material  instability.  When  Eulerian 
kernel  is  used,  the  material  instability  is  distorted  severely  as  shown  in  Fig.  1 .2(b). 

.As  a  compressible,  homogeneous  and  isotropic,  elastic  material  proposed  by  Blatz  and  Ko,  Blatz-Ko 
material  has  the  following  potential  function 

w{i„i„i,)=£(i2j-'+zj-s)  (1.2) 

where  J  =  and  //  is  a  constant. 
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Fig.  1.3  :  Stable  domain  from  theoretical  solution  compared  with  the  solution  from  PDE 


Fig.  1.4  :  Stable  domain  calculated  from  particle  methods 

Knowles  and  Sternberg  gave  the  stable  criteria  for  Blatz-Ko  material  by  an  analysis  of  the  associated 
issue  of  ellipticity.  The  stable  domain  is  a  hexagon  which  is  shown  in  Fig.  1 .4.  The  linearized  stability 
analysis  is  used  as  before  in  continuum  and  particle  methods  to  study  the  stable  domain,  compared  to  the 
theoretical  solution.  When  we  use  the  two-dimensional  stability  analysis  for  the  continuum,  the  stable 
domain  is  located  between  the  two  rays  from  the  original  point  as  shown  in  Fig.  1 .4.  It  is  identical  to  the 
stable  domain  in  Fig.  1 .4,  which  is  calculated  from  the  particle  methods  with  Lagrangian  kernel.  Fig.  1 .4 
also  shows  that  the  stable  domain  is  extremely  distorted  when  calculated  from  the  particle  methods  with 
Eulerian  kernel. 

2.  Effect  of  Cohesive  Laws  on  Crack' 

In  order  to  study  cohesive  cracks,  we  have  developed  a  new  enrichment  techniques  by  which  curved 
cohesive  cracks  can  be  treated  with  higher  order  enrichments.  The  crack  tip  can  be  located  anywhere 
within  an  element.  The  technique  is  based  on  Chen’s  (2003)  development  for  constant  strain  triangles 
and  extended  to  quadratic  elements.  An  essential  feature  in  applying  the  method  to  the  static  cohesive 
crack  model  is  that  the  crack  opening  should  be  quadratic  in  the  vicinity  of  the  crack  tip. 

Moes  and  Belytschko  (2002)  implemented  the  XFEM  scheme  for  cohesive  cracks  in  linear  triangular 
elements.  In  their  work,  a  linear  element  is  enriched  by  the  step  function  if  it  is  completely  cut  by  a  crack, 
and  by  the  branch  function  if  the  crack  tip  is  located  inside  the  element. 
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In  the  enrichment  developed  here,  an  element  containing  a  crack  tip  is  enriched  by  the  step  function 
whose  gradient  is  normal  to  the  crack  and  a  polynomial  along  the  crack.  We  consider  polynomials  of  first 
and  second  order  but  the  method  is  applicable  to  any  order.  The  contribution  of  the  paper  lies  in  that  the 
branch  function  used  in  previous  XFEM  (Moes  and  Belytschko  2002)  formulations  for  cracks  is  dropped, 
so  that  the  partition  of  unity  holds  in  the  entire  enriched  subdomain.  It  is  shown  that  this  yields  good  stress 
results  near  the  crack  tip  for  cohesive  crack  models.  Both  a  first  order  scheme  and  a  second  order 
scheme  are  developed  by  this  procedure.  The  scheme  can  be  applied  to  various  nonlinear  softening 
laws. 

Branch  functions,  typically,  rm  sin(«9/  2)  where  m  =  0.5,  2.0,  etc.,  have  been  used  for  the  enrichment  of 
the  tip  element  (Belytschko  and  Black  1 999;  Moes  and  Belytschko  2002;  Stazi  et  al.  2003).  For  a 
discussion  of  the  branch  function,  see  Belytschko  et  al.  (2001).  When  branch  functions  are  used  in 
conjunction  with  step  functions,  as  in  Moes  et  al.  (1 999),  the  partition  of  unity  property  does  not  hold  in  the 
elements  surrounding  the  tip  element.  The  enrichment  of  those  elements  is  a  local  partition  of  unity  and  it 
must  be  blended  to  the  rest  of  the  domain  for  optimal  performance  because  the  branch  function  does  not 
vanish  at  the  edges  of  the  tip  element.  Note  that  the  branch  function  is  not  a  piecewise  constant  function 
like  the  sign  function.  For  a  discussion  of  local  partitions  of  unity  and  blending,  see  Chessa  et  al.  (2003). 

In  order  to  avoid  the  difficulties  associated  with  the  branch  functions  in  the  tip  element,  a  new  enrichment 
has  been  developed.  This  procedure  is  similar  to  that  given  by  Chen  (2003)  but  it  is  generalized  to 
quadratic  elements  in  this  work.  It  can  be  applied  to  elements  of  any  order  without  difficulty. 


Fig.  2.1  :  A  two  dimensional  domain  containing  a  cohesive  crack  inside  the  domain. 

Consider  a  crack  Tc  in  a  two  dimensional  body  Q,  with  boundary  T  as  shown  in  Fig.  2.2.  We  denote  by 
x  the  spatial  coordinates.  The  crack  is  described  implicitly  by  the  level  set  function  <f>(x)  =  0  and  its  two 
end  points  Xj ,  I  =  1,2 .  Depending  on  the  cohesive  law,  a  traction  tc  is  applied  across  the  crack  surface 
near  the  tip. 

The  displacement  field  u  of  the  body  can  be  additively  decomposed  into  a  continuous  part  ucom  and  a 
discontinuous  part  udisc :  u(x)  =  ucml (x)  +  udisc(x) .  The  discontinuous  part  of  the  displacement  field  is 
approximated  by  standard  shape  functions  Nj  (x) ,  so 

ucoAx)=  IX  (2.1 ) 

l*Nm 

where  Nlol  is  the  total  set  of  nodes  and  u,  are  nodal  displacements. 
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Fig.  2.2  :  The  enriched  domain  Qenr  (shaded  area)  and  the  enriched  nodes  Nenr  (circles),  where  a 
virtual  (dashed)  line  segment  is  added  to  calculate  the  signed  distance  function  in  the  tip  element;  at 
square  nodes,  the  enrichment  parameters  a}  =  0 

If  we  consider  a  triangulation  on  the  body  Q.  as  shown  in  Fig.  2.2,  then  the  discontinuous  part  of  the 
displacement  field  can  be  limited  to  the  elements  that  contain  the  crack.  We  call  this  subdomain  Qenr . 

Let  Nmr  be  the  set  of  nodes  of  the  elements  cut  by  the  crack  rc .  The  discontinuous  part  of  the 
displacement  field  can  then  be  written  as  (Belytschko  and  Black  1999). 

«*,(*)=  'Zwm**,  (2-2) 

where  ^(x)  =  sign(<f>{x))  -  sign(0j )  and  a,  are  enrichment  parameters. 


Fig.  2.3  :  Enrichment  displacement  field  udisc  in  a  6-node  quadratic  tip  element  for  two  crack  directions 
and  the  corresponding  parent  domains;  the  shaded  areas  represent  the  enriched  region  Qenr . 

The  enrichment  developed  for  a  3-node  element  by  Chen  (2003)  is  generalized  to  a  6-node  element  (Fig. 
3)  as  follows.  We  consider  an  element  where  the  crack  passes  through  side  24.  Let  the  direction  of  the 
crack  be  such  that  it  intersects  side  13.  Other  relationships  between  the  element  and  the  crack  can  be 
obtained  by  permuting  the  node  numbers.  We  wish  to  construct  an  enrichment  that  vanishes  on  the  two 
edges  12  and  13  (which  is  necessary  for  compatibility),  and  is  continuous  across  edge  23  with  the  field  in 
the  adjacent  element.  To  meet  these  conditions,  only  nodes  3,  5,  6  are  enriched  and  the  discontinuous 
displacement  field  in  the  tip  element  is 
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*W  =  +  N5(f*)%(t*)a5  +  N6tf*)%tf*)a6  (2.3) 

where  £  are  the  parent  coordinates  of  the  shaded  region  23P  in  Fig.  3b,  i.e.  the  element  generated  by 
side  23,  side  3  P  and  P2,  where  P  is  the  intersection  point  of  the  line  joining  node  2  to  the  crack  tip  and 
side  32.  The  shaded  parent  area  coordinates  are  related  by  <f3’  =  1  -  -  £2*  and 

^  (£* )  =  sign{(f>{<!;‘ ))  -  sign(03 ) .  The  relation  between  t;  and  £  is  given  by 

C~4-  and  £=£  (2.4) 

fir 

where  %XP  is  the  area  coordinate  of  point  P. 

When  the  direction  of  the  crack  intersects  side  31  (Fig.  3c, d),  then  the  discontinuous  part  of  the 
displacement  field  is 

■W  =  +  W.  (£*)*,(?><  +  *,(#*)¥,(£*>,,  (2.5) 

where 

and  £=-&-  (2.6) 

b2P  b2  P 

Based  on  the  preceding,  we  can  simply  write  the  discontinuous  displacement  of  a  quadratic  tip  element  as 

•w-  (£>,(£)«,  m 

with  a,  =  0  at  three  of  the  nodes  as  shown  in  Fig.  2.3. 


Fig.  2.4  :  A  double  cantilever-beam  (DCB)  with  30%  notch  and  linear  and  bilinear  softening  laws. 

In  order  to  study  the  new  XFEM  scheme,  we  solved  the  double-cantilever-beam  (DCB)  problem  shown  in 
Fig.  2.4;  a  linear  softening  law  is  used.  Young’s  modulus  is  36.5  GPa  and  Poisson’s  ratio  is  0.18.  The 
results  obtained  by  the  new  XFEM  scheme  are  compared  to  those  obtained  by  the  static  condensation 
method  (e.g.  Zi  and  Bazant  2003)  which  has  extensively  been  used  to  analyze  the  fracture  properties  of 
many  experimental  tests  when  the  direction  of  crack  propagation  is  known.  The  static  condensation 
method  is  also  known  as  the  pseudo  boundary  integral  (PBI)  method. 

The  weak  form  and  the  discretized  equation  of  equilibrium  are  similar  to  Moes  and  Belytschko’s  (2002), 
so  it  is  omitted  here.  But,  the  condition  that  the  stress  at  the  crack  tip  should  be  equal  to  the  tensile 
strength  is  used  instead  of  the  zero  stress  intensity  factor  for  its  simplicity.  The  equilibrium  condition  and 
the  stress  condition  are  solved  simultaneously  by  the  Newton-Raphson  method. 
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Fig.  2.5  :  Relative  errors  in  load-deflection  histories  obtained  by  the  extended  finite  element  methods 
with  linear  elements  (empty  circles)  and  quadratic  elements  (solid  circles). 

The  convergence  rates  for  the  linear  and  quadratic  elements  for  various  mesh  refinements  are  shown  in 
Fig.  2.5  (left).  The  convergence  rate  of  the  second  order  method,  2.23,  is  much  better  than  the  first  order 
method,  0.48,  respectively. 

The  load-deflection  histories  for  the  quadratic  elements  with  the  stress  condition  are  plotted  and  compared 
to  the  result  by  the  static  condensation  method  in  Fig.  2.5  (right).  The  load  deflection  histories  are 
calculated  for  three  different  mesh  refinements  and,  for  all  cases,  agree  well  with  the  reference. 


line:  h=2  mm,  dash-dot  line:  h=5  mm,  dotted  line:  h=11  mm)  and  the  static  condensation  method 
(solid  line). 


The  crack  opening  and  the  stress  cry  is  shown  in  Fig.  2.6  for  three  different  mesh  refinements.  Because 

the  enrichment  function  is  based  on  a  partition  of  unity,  spurious  crack  opening  is  not  observed  even  in  the 
coarse  mesh.  The  crack  opening  profiles  for  h=  2  mm  almost  coincide  with  the  reference  crack  opening 
profiles  so  that  the  difference  is  indiscernible. 

The  stresses  on  the  dashed  line  in  Fig.  2.4(a)  are  plotted  in  Fig.  2.6.  It  can  be  seen  that  the  change  of  the 
stress  profiles  as  the  crack  advances  in  the  fracture  process  zone  is  approximately  self-similar.  As  the 
mesh  is  refined,  the  stress  profiles  converge  to  the  reference  profiles.  The  stress  profiles  for  h=2  mm 
agree  well  with  the  reference  profiles. 
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Fig.  2.7:  Comparison  of  the  results  by  the  linear  (solid  line)  and  the  bilinear  (dashed  line)  softening 
law;  the  load-deflection  histories  (left)  and  the  stresses  cry  (right). 

Because  of  their  simplicity,  linear  softening  laws  are  frequently  used.  However,  when  post-peak  behavior 
is  important,  nonlinear  softening  laws  may  be  more  accurate.  The  method  developed  can  be  applied  to 
nonlinear  softening  laws  as  well  as  linear  softening  laws.  One  might  think  that  if  the  fracture  energy  of  two 
different  softening  law  is  identical,  the  load-deflection  curve  does  not  change. 

In  order  to  study  the  effect  of  cohesive  laws,  a  bilinear  softening  law,  shown  in  Fig.  2.4  (right),  is  taken  into 
account.  The  tensile  strength  and  the  fracture  energy  of  the  bilinear  softening  law  (dashed  line)  are  the 
same  as  those  of  the  linear  softening  law  (solid  line).  The  load-deflection  histories  for  the  two  are 
compared  in  Fig.  2.7  (left).  The  peak  load  associated  with  the  bilinear  softening  law  is  smaller  than  the 
linear  one  and  the  shapes  of  the  load-deflection  histories  differ  as  well.  As  one  can  see  in  Fig.  2.7  (right), 
the  fracture  process  zones  are  not  fully  developed  at  the  peak  loads,  to  which  only  the  initial  parts  of  the 
softening  curves  in  Fig.  2.4  (right)  contribute.  Therefore  the  fracture  energies  associated  with  the  peak 
loads  are  different  although  the  fracture  energies  for  the  two  softening  laws  are  the  same. 


3.  Adaptivity  in  the  element  free  Galerkin  (EFG)  method 

We  have  studied  the  adaptive  procedures  that  provide  accurate  solutions  to  dynamic  nonlinear  with 
cohesive  cracking  problems  with  minimal  user  intervention.  The  key  feature  of  the  proposed 
methodologies  are  the  generation  of  models  with  a  prescribed  level  of  accuracy  via  adaptivity. 

The  raison  d’etre  for  meshfree  methods  is  to  totally  avoid  any  mesh  generation.  Thus  a  true  meshfree 
method  should  be  able  to  start  with  a  solid  model  or  geometric  model  and  provide  a  solution  of  the 
requisite  accuracy.  Obviously,  adaptivity  is  a  key  ingredient  in  any  such  procedure,  since  in  automatic 
model  generation,  user  intervention  in  the  assignment  of  refinement  would  be  unacceptable. 

The  procedure  proposed  here  combines  meshfree  methods  with  so-called  Cartesian  (or  structured 
methods).  In  these  methods,  the  arrangement  of  the  nodes  (or  particles)  is  structured  as  shown  in  Fig. 
3.2.  Although  at  first  such  a  combination  appears  illogical,  since  a  key  advantage  of  meshfree  methods  is 
their  ability  to  handle  arbitrary  clouds  of  nodes,  after  careful  consideration  we  concluded  there  is  little 
reason  to  model  objects  with  a  random  cloud  of  particles,  and  in  fact,  the  construction  of  random  clouds 
with  effective  fidelity  would  be  more  difficult  than  the  construction  of  structured  node  sets.  Furthermore, 
random  node  arrangements  would  complicate  the  implementation  of  adaptivity,  post-processing  and 
parallelization. 

One  could  question  the  use  of  meshfree  methods  in  a  structured  context;  if  the  nodal  arrangement  is 
structured,  it  would  appear  that  standard  finite  difference  methods  would  be  optimal.  However,  finite 
difference  methods  are  quite  unwieldy  in  treating  the  irregular  arrangements  of  nodes  that  arise  in  an 
adaptive  procedure.  Furthermore,  meshfree  methods  through  partition  of  unity  concepts  enable  the 
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method  to  incorporate  features  such  as  internal  discontinuities,  singular  solutions,  and  handbook 
solutions.  These  features  endow  the  methodology  with  powerful  advantages:  discontinuous  partitions  of 
unity  can  treat  arbitrary  cracks  without  remeshing,  singular  partitions  of  unity  can  treat  elastic  near  tip 
fields  for  cracks  and  dislocations  without  mesh  refinement,  handbook  solutions  can  treat  holes  and 
inclusions  that  are  smaller  than  the  scale  of  the  discretization. 

Some  of  these  ideas  have  already  been  presented  in  the  finite  element  context.  Strouboulis  et  al.  (2001) 
in  the  generalized  finite  element  method  introduced  the  concept  of  adding  handbook  solutions  for 
subscale  features.  Belytschko  et  al.  (2003)  developed  a  structured  finite  element  method  based  on  implicit 
surface  definitions  with  arbitrary  internal  features  such  as  cracks  and  and  sliding  interfaces.  Belytschko  et 
al.  (2003)  showed  how  this  methodology  could  be  applied  to  topology  optimization. 

Another  question  the  reader  may  pose  is  why  use  meshfree  methods  instead  of  higher  order  finite 
elements.  Indeed  spectral  finite  elements  are  extremely  powerful  for  treating  many  problems.  However,  for 
the  nonsmooth  problems  found  in  solid  mechanics,  such  aselastic-plastic  materials  and  contact-impact,  p- 
adaptivity  is  inferior  to  h-adaptivity.  Furthermore,  the  implementation  of  h-adaptivity  for  finite  elements 
even  in  a  structured  mesh  is  much  more  difficult  than  for  meshfree  methods. 

We  will  use  the  element-free  Galerkin  (EFG)  meshfree  method  in  this  study.  This  method  uses  moving 
least  square  (MLS)  approximations.  We  will  not  consider  some  recent  methods  such  as  radial  basis 
functions  and  so  called  truly  meshless  methods.  The  former  require  multipolar  methods  for  reasonable 
speed  on  large  models.  The  truly  meshless  methods  are  more  in  the  spirit  of  clouds  of  particles,  which  are 
not  in  the  spirit  of  the  methodology  proposed  here. 

We  have  chosen  a  Total  Lagrangian  description.  The  momentum  equation  can  then  be  written  as 

\paSu,b,da. o-  J&^dTo-O  (3.1) 

n0  o0  dX‘  flo  n« 

where  /?„  is  the  initial  density,  t/(  are  the  displacements,  P*  are  the  nominal  stresses,  b,  are  the  body 
forces  and  l  are  the  external  forces.  The  test  functions  are  approximated  by 

Sut^&taOjiX,)  (3.2) 

JzS 


where  the  shape  functions  <bj(Xt)  are  the  EFG  shape  functions.  The  integrals  are  evaluated  by  an 
integration  based  on  a  background  mesh. 

The  adaptivity  approach  is  based  on  an  interpolation  error.  Regions  with  a  large  local  error,  i.e.  high  strain 
gradients,  are  refined.  The  refinement  strategy  is  illustrated  in  Fig.  3.2.  The  cell  is  divided  into  four  smaller 
cells  of  equal  size.  In  our  numerical  studies  large  errors  occurred  if  the  size  of  the  neighboring  cells  differs 
too  much.  Therefore,  a  cell  is  refined,  if  it  shares  more  than  two  edges  with  a  neighbor  cell  as  shown  in 
Fig.  3.2  where  the  last  refinement  step  from  Fig.  3.1  is  chosen  as  source. 

In  dynamic  problems,  a  lumped  or  diagonal  mass  matrix  is  used  instead  of  a  consistent  one  since  it 
simplifies  the  computation.  In  particle  methods,  there  exist  different  techniques  to  compute  the  masses  of 
the  particles.  A  common  procedure  is  to  assign  the  masses  and  volumes  to  the  particles  according  to  a 
Voronoi  diagram  (see  Fig.  3.3)  before  beginning  the  time  integration.  Another  way  is  to  establish  a 
diagonal  mass  matrix  from  the  consistent  mass  matrix  using  the  row  sum  technique. 

If  adaptivity  takes  place,  it  becomes  necessary  to  recompute  the  new  masses  and  volumes  after  every 
refinement  step.  It  is  computationally  easier  to  compute  the  diagonal  mass  matrix  from  the  consistent  one, 
especially  if  particles  are  added  in  an  irregular  pattern.  Therefore,  we  have  chosen  this  approach. 
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If  an  incremental  stress  strain  law  is  applied,  the  stresses  have  to  be  assigned  to  the  new  particles.  Also 
values  of  the  velocity  have  to  be  transformed  to  the  new  particles.  All  values  are  obtained  by  an  MLS  fit, 
i.e.  via  the  approximation.  Especially,  when  the  particles  are  arranged  in  an  irregular  pattern,  this 
procedure  is  straightforward.  The  new  interpolation  radius  h  is  chosen  according  to  the  maximum  distance 
of  the  neighbor  particles.  In  the  particle  sums,  an  average  of  hu  =0.5(/z;  +hj) is  calculated  to  obtain  a 

symmetric  interaction. 

The  response  of  notched  concrete  under  dynamic  loading  have  been  studied.  Here,  we  consider  the  beam 
as  shown  in  Fig.  3.4.  We  compare  again  the  crack  pattern  obtained  with  a  continuum  model  and  the 
discrete  crack  model.  EFG  with  cell  integration  was  taken  to  solve  this  problem.  Fig.  3.5  shows  the  crack 
obtained  from  the  simulation  with  a  continuum  damage  model  and  the  discrete  crack  model.  It  can  be 
seen  that  the  adaptive  method  refines  around  the  path  of  the  crack,  so  that  a  very  accurate  path  can  be 
obtained.  Fig.  3.  6  illustrated  the  refinement  of  the  crack  close  to  the  notch. 


Fig.  3.2  :  Refinement  to  create  a  smooth  transition 


•  Miaster  particles 


Fig.  3.3  :  Voronoi  diagram 
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Fig.  3.4  :  The  notched  concrete  beam  of  John  and  Shah 


Fig.  3.5  :  The  crack  pattern  of  the  notched  concrete  beam  in  the  numerical  simulation 
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Fig.  3.6  :  Refinement  around  the  crack 


4.  Domain  decomposition 


A  method  has  been  developed  for  coupling  fine  scale  models  of  subdomains  with  coarse  scale  models. 
The  essential  feature  of  any  such  coupling  is  that  it  should  not  lead  to  spurious  reflections  off  the  interface. 
Direct  linkages  of  fine  and  coarse  scale  models  lead  to  nonphysical  reflections  because  the  impedance  of 
the  coarse  scale  model  differs  significantly  from  that  of  the  fine  scale  model.  In  face,  for  frequencies 
above  the  cutoff  frequency  of  the  fine  scale  model,  the  coarse  scale  model  appears  almost  as  a  rigid 
boundary  condition. 


The  coupling  method  is  based  on  overlapping  domain  decomposition  techniques,  see  Belytschko  and 
Xiao  (2003).  The  essential  idea  is  that  on  the  overlapping  domain,  the  fine  scale  solution  is  projected  in 
each  time  step  onto  the  coarse  scale  solution..  This  is  accomplished  by  scaling  the  Hamiltonian  within  the 
bridging  domain  and  constraining  the  velocities  of  the  overlapping  nodes  by  Lagrange  multipliers. 

The  Hamiltonian  is  defined  as: 


r  P2 

H  =  K+E  =  L  — 
2  P 


ii 


W(F)dCl 


(4.1) 


where  P  is  the  momentum,  p  is  the  density,  F  is  the  gradient  of  deformation,  and  W  is  the  internal 
energy  or  strain  energy. 


The  equations  of  motion  can  be  written  as: 


and  P 


m 

dx 


(4.2) 


The  total  Hamiltonian  of  Bridging  domain  multiscale  method  for  FEM  is  written  as  the  linear  combination  of 
Hamiltonian  from  fine  mesh  QF  and  coarse  mesh  Qc  of  FEM,  as  shown  in  Fig.  4.2. 

H  =  pHF  +  (1  -  p)Hc  =  p(KF+  Ef  )+  (l  -  piKc  +  Ec )  (4.3) 

where  P  is  a  scalar  parameter  which  ranges  from  0  to  1  in  bridging  domain  Qint .  The  bridging  domain  is 
also  called  overlapping  domain. 
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Fig.  4.1  :  Overlapping  domain  of  bridging  domain  muttiscale  method 


The  two  models  are  constrained  on  the  overlapping  subdomain  Qm<  by 

(X,)-rf,,}  =  few 

i.e.  the  fine  mesh  displacements  d  conform  to  the  coarse  mesh  displacements  u  at  the  discrete 
positions  of  the  fine  mesh  nodes.  The  constraints  are  applied  on  the  each  component  of  the 
displacements. 

When  Lagrangian  multiplier  method  is  applied,  the  total  Hamiltonian  will  be  given  as: 

HL=H  +  XTg  =  H  +  YiX]  g,  (4.5) 

i 

where  X  is  the  vector  of  Lagrange  multipliers  which  are  assigned  on  the  discrete  fine  mesh  nodes  in  the 
bridging  domain. 


g  i  =  fc/}=  k 


.Or,)*-*-. 


(4.4) 


Then,  the  equations  of  motion  are: 

p  =  _?!h 
dx 
i.e. 

rrr  _  i -ext  jMnt  Y)  T  _  *  ext  i*int  *L  .  nf 

m/d;  -  fj  -f;  '"Q 

M,ii,  =  F"  -  F;nl  -  Y  Xj  —  =  Ff  -  Fj“  -  Fj-  in  Clc 

j  ^U/ 

where 

(x)w,(x)«0F 

m,  =  E  0  -  Mx)]N,{x^j{x)dac0 

«r  - 1 

F"  =  ifO-Ax))^^p<n? 


»/  =Ei fAx)w 


(4.6) 
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tf  =  fi(X}N,p0b^  +  [Ax)N,ldTl 

F"  =  +  ^(l-^(X)Kt< 


We  developed  an  explicit  algorithm  based  on  central  difference  method  to  solve  these  equations  of  motion. 
Suppose  we  already  know  the  accelerations,  velocities  and  displacement  for  time  step  n .  We  first  obtain 
the  displacements  at  the  next  time  step  n  + 1 . 


U/(n+l) 

•  *  1 

=  um+nmAt  +  -i 

in 

a 

^/(n+1) 

=  d/(«)  +  <*/(„)  Af  +  2* 

^/(»)A?2 

in 

Q 

(4.7) 


M 

0 


The  accelerations  for  next  time  step  can  be  given  from  (4.6)  without  considering  the  force  due  to  the 
constraints: 

1 


*/(«+]) 


u 


/(»+!) 


m, 


M, 


|  next  fint 

^/(n+l)  “I/(w+l) 


] 


-r?ex/  pi  in/  1 

*/(/ 7+1)  “r/(*+l)J 


in  Q ; 


in  Qq 


(4.8) 


Then,  we  obtain  the  trial  velocities: 


ll(n+ 1)  -  u/(n) 


■3 


+  —  d  rr„\  +  d 


/(it)  ^  w/(w+l) 


”/(»+!)  =  */(»)  +^t»/(n)  +  »/(«+!) 


in  Q 


0 


w  Qq 


(4.9) 


Alternatively,  the  velocities  at  time  step  «  +  1  can  be  expressed  as: 


J/(n+l)  -  d/(»)  +  —  [d/(n)  -  m}  +  d7(n+1)  -  mj  !/(„+!) 


--UL 


— lfi 

-  m,  f, 


h 


u 


-  <*/(„+!)  -  mj  'A /£g5*, 

J 

/(«+ 1)  =  “/(„)  +  ^  [**/(«)  -  Mj  F/w  +  U/(n+1) 


-  “/(n+l)  ~  Ml 


where  G^j 


dg/ 
du , 


[wJti] ,  g£  = 


dgl 

dd7 


(4.10) 


(4.11) 


-  [“  ^/zl]  ■  *-/  -  T  [^7(«)  +  ^/(«+i)  ]  denote  the 


unknown  Lagrange  multipliers.  The  above  velocities  satisfy  the  constraints  (4)  in  time  derivative  form 

8 /(»+!)  =  “(X/ )(„+i)  1)  =  X^/(Xj)“y(«+l)  -^7(«+l)  (4.12) 


Substituting  (4.10)  and  (4.11)  into  (4.12),  the  unknown  Lagrange  multipliers  can  be  obtained  by  solving  the 
following  equations: 

-  8/  (4.13) 

L 
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where 


A/£  -  AtMj  l^dNjIG^j  ,  g }  -  y^Njjiij  d7 

j  j 

After  obtaining  the  Lagrange  multipliers,  the  corrected  velocities  at  time  n  +  l  are  given  in  (4.10)  and 
(4.11). 

To  study  the  effectiveness  of  this  method,  a  wave  propagation  problem  in  a  one-dimensional  rod  is 
studied  for  both  the  bridging  domain  coupling  method  and  an  edge-to-edge  coupling  method.  The  initial 
wave  is  a  combination  of  high  frequency  wave  and  low  frequency  wave.  Spurious  wave  reflection  occurs 
when  edge-to-edge  coupling  method  is  used  as  shown  in  Fig.  4.2(b).  However,  Fig.  4.2(a)  shows  that 
bridging  domain  coupling  method  can  eliminate  the  wave  reflection  without  any  additional  filtering  or 
viscosities.  Fig.  4.3  shows  the  histories  of  energy  in  fine  mesh  domain  or  coarse  mesh  domain.  After  the 
low  frequency  wave  passes  through  the  coarse  mesh  domain,  similar  energies  pass  through  whatever 
coupling  method  is  used.  In  fine  mesh  domain,  the  energy  of  high  frequency  remains  from  edge-to-edge 
coupling  method  due  to  the  spurious  wave  reflection,  but  it  is  reduced  to  zero  from  bridging  domain 
coupling  method  as  shown  in  Fig.  4.3(a). 


Fig.  4.2  :  Spurious  wave  reflection  in  multiscale  methods 


(b)  coarse  mesh  domain 
Fig.  4.3  :  Histories  of  energy 


Cylindrical  wave,  propagation  in  two-dimensional  plate  is  shown  in  Fig.  4.4.  Fig.  4.4(a)  shows  the  bridging 
domain  coupling  model  which  we  use  here.  Fig.  4.4  (b)-(d)  show  the  configurations  of  wave  propagation  at 
three  different  times.  There  is  almost  no  spurious  wave  reflection  when  cylinder  wave  propagates  from 
fine  mesh  domain  to  coarse  mesh  domain. 
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^flne  ^overlap  ^coarse 

(a)  bridging  domain  model  (b)  t  =  1  second 


(c)  t  =  2  second  (d)  t  =  3  second 

Fig.  4.4  :  Cylinder  wave  propagation  example 
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